Frecuentismo: ¿qué edad tenéis?

• La media de edad es una cantidad fija, pero desconocida (\(\mu_{\text{edad}}\)).
• El objetivo es estimar \(\mu_{\text{edad}}\) usando una muestra aleatoria.
  1. Elegimos un estimador: la media muestral.
  2. Cuantificamos la incertidumbre: el error típico de la media (SE).
  3. Usamos el SE y construimos un intervalo de confianza (IC) para inferir el valor de \(\mu_{\text{edad}}\).
• Un frecuentista concluirá:
  • La media de edad estimada es \(\hat\mu_{\text{edad}}\) con un IC de \([a,b]\).
  • Si repetimos este experimento infinitas veces, el IC contendrá el valor de \(\mu_{\text{edad}}\) en el 95% de los casos.

Bayesiano: ¿qué edad tenéis?

• La media de edad es una variable aleatoria porque incluimos nuestra incertidumbre sobre su valor con una distribución de probabilidad.
• El objetivo es combinar (1) los datos y (2) lo que creemos.
  1. Elegimos una distribución a priori para \(\mu_{\text{edad}}\): \(\Pr(\mu_{\text{edad}})\)

Bayesiano: ¿qué edad tenéis?

Bayesiano: ¿qué edad tenéis?

• Un bayesiano concluirá:
  • Tras ver los datos, mi antigua creencia (prior) sobre \(\mu_{\text{edad}}\) se ha actualizado (posterior) .
  • El intervalo de credibilidad al 95% nos dice que hay un 95% de probabilidad de que \(\mu_{\text{edad}}\) esté en \([a,b]\)

Universo frecuentista, multiverso bayesiano

Universo frecuentista, multiverso bayesiano

Universo frecuentista, multiverso bayesiano

Universo frecuentista, multiverso bayesiano

Cada año, Carmen Ximénez y Javier Revuelta, como coordinadores del Máster de Metodología de las Ciencias del Comportamiento y de la Salud, deben decidir cuántos Trabajos de Fin de Máster (TFM) de investigación y de prácticas van a ofertar.

Para ello, necesitan conocer la proporción esperada de estudiantes que optarán por un TFM de investigación. Denotemos esta proporción como \(\theta\), un valor desconocido que podemos estimar a partir de los datos del curso anterior.

# TFMs de investigación el último año
TFM <- c(1,0,0,1,1,1,0,1,1,0,
         0,0,1,0,1,1,1,1,0,1)

Estimando \(\theta\): método de los momentos

En esta situación, el método de los momentos es la media aritmética de la variable TFM.

# Estimación de \theta con el método de los momentos
mean(TFM)

[1] 0.6

Como \(\hat\theta=0.6\), el 60% de los alumnos escogieron un TFM de investigación el año pasado.

Estimando \(\theta\): máxima verosimilitud

• Necesitamos un modelo de probabilidad adecuado.
  • Dado que cada estudiante elige entre dos opciones, podemos modelar su decisión con la distribución de Bernoulli

• Probabilidad de escoger un TFM de investigación: \(\theta\)
• Probabilidad de escoger un TFM de prácticas: \(1-\theta\)
• En R, la ecuación anterior se representa con dbinom, indicando size=1.

# Función de densidad de bernoulli
dbinom(x = "TFM_i", size = 1, prob = "theta")

Estimando \(\theta\): máxima verosimilitud

• Sin embargo, nuestro interés no está en cada estudiante por separado, sino en la verosimilitud conjunta de los datos.
• Cada alumno elige un TFM independientemente de lo que elijan los demás.
• La verosimilitud conjunta se obtiene multiplicando las verosimilitudes individuales de los alumnos.

# Función de verosimilitud de los datos
prod(dbinom(x = "TFM_i", size = 1, prob = "theta"))

Estimando \(\theta\): máxima verosimilitud

• ¿Cómo de verosímil es que \(\theta=0.2\) en nuestros datos?

# Verosimilitud de los datos cuando \theta = 0.2
prod(dbinom(TFM, size = 1, prob = 0.2))

[1] 6.871948e-10

• El valor 6.871948e-10 representa la verosimilitud de los datos dado un valor de \(\theta = 0.2\). ¿Qué tal si probamos ahora \(\theta = 0.3\)?

# Verosimilitud de los datos cuando \theta = 0.3
prod(dbinom(TFM, size = 1, prob = 0.3))

[1] 3.063652e-08

• Como 3.063652e-08 es mayor que 6.871948e-10, \(\theta=0.3\) explica mejor los datos que \(\theta=0.2\).

Estimando \(\theta\): máxima verosimilitud

¿Qué pasa si repetimos este mismo proceso para todos los valores de \(\theta\) entre 0 y 1?

# Valores de theta de 0 a 1
theta <- seq(0, 1, 0.01)
# Vector vacío para guardar la verosimilitud con cada theta
likelihood <- vector(length = length(theta))
# Verosimilitud para cada theta
for(i in 1:length(theta)){
  likelihood[i] <- prod(dbinom(TFM, size = 1, prob = theta[i]))
}
# Gráfico de densidad
plot(theta, likelihood)

Estimando \(\theta\): máxima verosimilitud

• ¿Cuál es el valor de \(\boldsymbol\theta\) que maximiza la verosimilitud?

Estimando \(\theta\): máxima verosimilitud

Estimando \(\theta\): máxima verosimilitud

Estimando \(\theta\): máxima verosimilitud

# Valor de theta que maximiza la verosimilitud
theta[which.max(likelihood)]

[1] 0.6

• Este valor de \(\hat\theta=0.6\) coincide con el que hemos obtenido antes utilizando el método de los momentos.

Estimando \(\theta\): estimación bayesiana

• La estimación bayesiana es, ante todo, estimación basada en la verosimilitud, y su corazón es el Teorema de Bayes:

• Verosimilitud: \(\Pr\left(D \mid \theta\right)\)
• Probabilidad a priori: \(\Pr(\theta)\)
• Probabilidad posterior: \(\Pr\left(\theta \mid D\right)\)

Estimando \(\theta\): estimación bayesiana

• Alternativamente, el Teorema de Bayes también se representa como

• \(\Pr(D)\) recibe el nombre de verosimilitud marginal, y no es más que el sumatorio de todos los valores de \(\Pr\left(D \mid \theta\right) \times \Pr(\theta)\).

# Ejemplo
likelihood <- c(.2, .8, .8); prior <- c(.2, .5, .5)
# Sin dividir entre su sumatorio
likelihood * prior

[1] 0.04 0.40 0.40

Estimando \(\theta\): estimación bayesiana

# Ejemplo
likelihood <- c(.2, .8, .8); prior <- c(.2, .5, .5)
# Dividiendo entre su sumatorio
likelihood * prior/sum(likelihood * prior)

[1] 0.04761905 0.47619048 0.47619048

Estimando \(\theta\): estimación bayesiana

• Esto garantiza que la distribución posterior, \(\Pr\left(\theta \mid D\right)\), sea una distribución de probabilidad válida.
• La verosimilitud marginal es imprescindible en la comparación de modelos a través del factor de Bayes.

Estimando \(\theta\): estimación bayesiana

• Carmen y Javier han recuperado la selección de TFMs de los cursos académicos.

• \(\hat\theta\) es nuestro conocimiento previo de la selección de TFMs.

Estimando \(\theta\): estimación bayesiana

• La media de \(\hat\theta\) es 0.40, y su desviación típica es 0.12.
• Podemos utilizar esta información con la distribución a priori.

Estimando \(\theta\): estimación bayesiana

• El siguiente código de R muestra cómo surge toda la magia que crea la famosa distribución posterior.
• En primer lugar, calculamos (otra vez) la verosimilitud de los datos

# Valores de theta de 0 a 1
theta <- seq(0, 1, 0.01)
# Vector vacío para guardar la verosimilitud con cada theta
likelihood <- vector(length = length(theta))
# Verosimilitud para cada theta
for(i in 1:length(theta)){
  likelihood[i] <- prod(dbinom(TFM, size = 1, prob = theta[i]))
}

Estimando \(\theta\): estimación bayesiana

• Aplicamos el Teorema de Bayes:

# Distribución a priori de theta
prior <- dnorm(theta, mean = 0.4, sd = 0.12)
# Distribución posterior (proporcional)
posterior_prop <- likelihood * prior
# Verosimilitud marginal
marginal <- sum(likelihood * prior)
# Distribución posterior escalada a [0,1]
posterior <- posterior_prop/marginal 
plot(theta, posterior) # Gráfico de densidad posterior

Estimando \(\theta\): estimación bayesiana

\[ \Pr\left(\theta \mid D\right) = \frac{\Pr\left(D \mid \theta\right) \times \color{darkorange}{\Pr(\theta)}}{\Pr(D)}=\frac{\text{Verosimilitud }\times\color{darkorange}{\text{ Prior}}}{\text{Verosimilitud marginal}} \]

# Distribución a priori de theta
prior <- dnorm(theta, mean = 0.4, sd = 0.12)
# Distribución posterior (proporcional)
posterior_prop <- likelihood * prior
# Verosimilitud marginal
marginal <- sum(likelihood * prior)
# Distribución posterior escalada a [0,1]
posterior <- posterior_prop/marginal 
plot(theta, posterior) # Gráfico de densidad posterior

Estimando \(\theta\): estimación bayesiana

\[ \Pr\left(\theta \mid D\right) = \frac{\color{darkorange}{\Pr\left(D \mid \theta\right) \times \Pr(\theta)}}{\Pr(D)}=\frac{\color{darkorange}{\text{Verosimilitud }\times\text{ Prior}}}{\text{Verosimilitud marginal}} \]

# Distribución a priori de theta
prior <- dnorm(theta, mean = 0.4, sd = 0.12)
# Distribución posterior (proporcional)
posterior_prop <- likelihood * prior
# Verosimilitud marginal
marginal <- sum(likelihood * prior)
# Distribución posterior escalada a [0,1]
posterior <- posterior_prop/marginal 
plot(theta, posterior) # Gráfico de densidad posterior

Estimando \(\theta\): estimación bayesiana

\[ \Pr\left(\theta \mid D\right) = \frac{\Pr\left(D \mid \theta\right) \times \Pr(\theta)}{\color{royalblue}{\Pr(D)}}=\frac{\text{Verosimilitud }\times\text{ Prior}}{\color{royalblue}{\text{Verosimilitud marginal}}} \]

# Distribución a priori de theta
prior <- dnorm(theta, mean = 0.4, sd = 0.12)
# Distribución posterior (proporcional)
posterior_prop <- likelihood * prior
# Verosimilitud marginal
marginal <- sum(likelihood * prior)
# Distribución posterior escalada a [0,1]
posterior <- posterior_prop/marginal 
plot(theta, posterior) # Gráfico de densidad posterior

Estimando \(\theta\): estimación bayesiana

\[ \color{red}{\Pr\left(\theta \mid D\right)} = \frac{\Pr\left(D \mid \theta\right) \times \Pr(\theta)}{\Pr(D)}=\frac{\text{Verosimilitud }\times\text{ Prior}}{\text{Verosimilitud marginal}} \]

# Distribución a priori de theta
prior <- dnorm(theta, mean = 0.4, sd = 0.12)
# Distribución posterior (proporcional)
posterior_prop <- likelihood * prior
# Verosimilitud marginal
marginal <- sum(likelihood * prior)
# Distribución posterior escalada a [0,1]
posterior <- posterior_prop/marginal 
plot(theta, posterior) # Gráfico de densidad posterior

Estimando \(\theta\): estimación bayesiana

Estimando \(\theta\): estimación bayesiana

• La moda posterior suele utilizarse como estadísticos de resumen de la distribución posterior.

# Moda de la distribución posterior
theta[which.max(posterior)]

[1] 0.51

• Es distinto del valor que obtendríamos por máxima verosimilitud

# Estimación máximo-verosimil
theta[which.max(likelihood)]

[1] 0.6

• La distribución a priori es responsable de esta diferencia

Estimando \(\theta\): estimación bayesiana

• ¿Qué pasa si utilizamos una distribución a priori menos informativa?

Estimando \(\theta\): estimación bayesiana

\[ \theta \sim \mathcal{N}\left(\mu_\theta = \color{royalblue}{0.5},\, \sigma_\theta = \color{royalblue}{0.1}\right) \]

Estimando \(\theta\): estimación bayesiana

\[ \theta \sim \mathcal{N}\left(\mu_\theta = \color{royalblue}{0.5},\, \sigma_\theta = \color{royalblue}{1}\right) \]

Estimando \(\theta\): estimación bayesiana

\[ \theta \sim \mathcal{N}\left(\mu_\theta = \color{royalblue}{0.5},\, \sigma_\theta = \color{royalblue}{5}\right) \]

Estimando \(\theta\): estimación bayesiana

\[ \theta \sim \mathcal{N}\left(\mu_\theta = \color{royalblue}{0.5},\, \sigma_\theta = \color{royalblue}{5}\right) \]

Estimando \(\theta\): estimación bayesiana

• La probabilidad a priori para cualquier valor de \(\theta\) es prácticamente la misma.

# Prior no informativa sobre theta
theta <- seq(0, 1, .1)
round(dnorm(theta, 0.5, 5), 4)

 [1] 0.0794 0.0795 0.0796 0.0797 0.0798 0.0798 0.0798 0.0797 0.0796 0.0795
[11] 0.0794

Estimando \(\theta\): estimación bayesiana

• Las distribuciones a priori planas, como la distribución uniforme, son los casos más extremos de distribuciones a priori no informativas.

• Se asigna una probabilidad constante a todo valor en el rango \([a,b]\)

Estimando \(\theta\): estimación bayesiana

theta <- seq(0, 1, .1) # Valores de theta
dunif(theta, min = 0, max = 1) # Prior uniforme

 [1] 1 1 1 1 1 1 1 1 1 1 1

• ¿Qué implicaciones tiene esto para el Teorema de Bayes?

Estimando \(\theta\): estimación bayesiana

theta <- seq(0, 1, .1) # Valores de theta
dunif(theta, min = 0, max = 1) # Prior uniforme

 [1] 1 1 1 1 1 1 1 1 1 1 1

\[ \Pr\left(\theta \mid D\right) \propto \Pr\left(D \mid \theta\right) \times \color{red}{\boldsymbol{1}}\\ \]

Estimando \(\theta\): estimación bayesiana

theta <- seq(0, 1, .1) # Valores de theta
dunif(theta, min = 0, max = 1) # Prior uniforme

 [1] 1 1 1 1 1 1 1 1 1 1 1

\[ \Pr\left(\theta \mid D\right) \color{red}= \Pr\left(D \mid \theta\right) \]

• La distribución posterior es lo mismo que la verosimilitud.

Máxima verosimilitud Vs. Cadenas de Markov

• Máxima verosimilitud utiliza algoritmos de optimización numérica para encontrar un valor máximo.

Máxima verosimilitud Vs. Cadenas de Markov

• En ocasiones, puede encontrar un máximo local, pero no un máximo absoluto.

Máxima verosimilitud Vs. Cadenas de Markov

• Las Cadenas de Markov utilizan métodos de simulación Monte Carlo para recuperar toda una distribución, no un único valor

Máxima verosimilitud Vs. Cadenas de Markov

Máxima verosimilitud Vs. Cadenas de Markov

Máxima verosimilitud Vs. Cadenas de Markov

• Los algoritmos que se suelen aplicar difieren en eficiencia y flexibilidad.
  • Sampleador de Gibbs (implementado en JAGS).
  • Monte Carlo Hamiltoniano (implementado en Stan).
• Visualmente podemos ver cómo difieren en su forma de extraer muestras.
  • Visualizador de algoritmos MCMC

¿Qué es modelar?

¿Qué es modelar?

¿Qué es modelar?

¿Qué es modelar?

• Es habitual expresar un modelo de regresión lineal como

• Esto es idéntico a

• Por tanto, el modelo estadístico cuenta con:
  • Una variable dependiente y una distribución de probabilidad,
  • los parámetros de la distribución de probabilidad seleccionada.
  • los pronósticos de la ecuación lineal (i.e., la parte determinista).

El modelo psicométrico del factor común

• La psicometría se asienta sobre el modelo de regresión lineal.

• Sin embargo, en psicometría:
  • se asume que la variable \(X\) es latente en lugar de observable. Se suele representar con la letra griega \(\eta\).
  • suelen utilizarse las letras griegas \(\nu\), \(\lambda\) y \(\theta\) en lugar de \(\beta_0\), \(\beta_1\) y \(\sigma_\varepsilon\).
• Por eso, el mismo modelo se representa como

El modelo psicométrico del factor común

\[ Y_i = \nu + \lambda \cdot \eta_i + \varepsilon_i,\quad\text{donde}\quad\varepsilon_i\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta\right) \]

• \(Y_i\) es una medida imperfecta de una variable no observable, \(\eta_i\).
• \(\varepsilon_i\) recibe el nombre de unicidad (i.e., varianza no explicada).
  • Una parte de la unicidad es el error de medida
• Es habitual contar con varias medidas de la variable no observable.

El modelo psicométrico del factor común

\[ Y_i = \nu + \lambda \cdot \eta_i + \varepsilon_i,\quad\text{donde}\quad\varepsilon_i\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta\right) \]

El modelo psicométrico del factor común

\[ Y_i = \nu + \lambda \cdot \eta_i + \varepsilon_i,\quad\text{donde}\quad\varepsilon_i\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta\right) \]

El modelo psicométrico del factor común

\[ Y_i = \nu + \lambda \cdot \eta_i + \varepsilon_i,\quad\text{donde}\quad\varepsilon_i\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta\right) \]

El modelo psicométrico del factor común

\[ Y_i = \nu + \lambda \cdot \eta_i + \varepsilon_i,\quad\text{donde}\quad\varepsilon_i\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta\right) \]

• ¿Y si tenemos más de una variable latente?

El modelo psicométrico del factor común

• ¿Cómo cambia el modelo al ir añadiendo variables latentes?

El modelo psicométrico del factor común

El modelo psicométrico del factor común

El modelo psicométrico del factor común

• Cuando hay \(M\) variables latentes, el modelo se puede escribir como

El modelo psicométrico del factor común

• Los psicómetras prefieren la versión matricial del modelo.

• En la versión matricial, para cada sujeto \(i\) tenemos:
  • Un vector con \(J\) respuestas observadas: \(\mathbf{y}_i\).
  • Un vector con \(J\) intersecciones: \(\boldsymbol\nu\).
  • Una matriz con \(J\times M\) pesos factoriales: \(\mathbf{\Lambda}\).
  • Un vector con \(M\) puntuaciones latentes: \(\boldsymbol\eta_i\).
  • Un vector con \(J\) puntuaciones error: \(\boldsymbol{\varepsilon}_i\).

El modelo psicométrico del factor común

• Ambas versiones representan dos formas distintas del mismo modelo
• La versión condicional

• La versión marginal

• \(\boldsymbol{\alpha}\) es el vector con las medias de \(\boldsymbol{\eta}_i\). \(\mathbf\Psi\) y \(\mathbf\Theta\) son las matrices de covarianzas de \(\boldsymbol{\eta}_i\) y \(\boldsymbol{\varepsilon}_i\), respectivamente.

El modelo psicométrico del factor común

\[ Y_{ij} \sim \mathcal{N}\left(\mu = \nu_j + \sum^M_{m=1}\lambda_{jm} \cdot \color{firebrick}{\eta_{im}},\,\sigma = \theta\right) \]

\[ \mathbf{y}_i \sim \mathcal{MVN}\left(\boldsymbol\mu = \boldsymbol\nu + \mathbf{\Lambda} \cdot \color{firebrick}{\boldsymbol{\alpha}},\,\mathbf{\Sigma} = \mathbf{\Lambda\cdot\color{firebrick}{\mathbf\Psi}\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}\right) \]

Simulando datos desde el modelo en R

Simulando datos desde el modelo en R

\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]

• Intersecciones de los ítems: \(\boldsymbol{\nu}\)

# Vector de intersecciones poblacional
Nu <- c(-1, 0, 1, -1, 0, 1)

Simulando datos desde el modelo en R

\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]

• Pesos factoriales no estandarizados: \(\mathbf{\Lambda}\)

Lambda <- matrix(c(1.40, 1.75, 2.10, 0, 0, 0,
                   0, 0, 0, 1.40, 1.75, 2.10), ncol = 2)
t(Lambda) # Lambda transpuesta 

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]  1.4 1.75  2.1  0.0 0.00  0.0
[2,]  0.0 0.00  0.0  1.4 1.75  2.1

Simulando datos desde el modelo en R

\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]

• Medias de las variables latentes: \(\boldsymbol{\alpha}\)

# Medias latentes poblacionales
alpha <- c(0, 0)

Simulando datos desde el modelo en R

\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]

• Matriz de covarianzas entre las variables latentes: \(\mathbf{\Psi}\)

# Matriz de covarianzas latente poblacional
Psi <- matrix(c(1, .5, .5, 1), ncol = 2)

Simulando datos desde el modelo en R

\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]

• Matriz de varianzas únicas: \(\mathbf{\Theta}\)

Simulando datos desde el modelo en R

\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]

\[ \mathbf{\Theta}=\begin{bmatrix} 2.04 & 0 & 0 & 0 & 0 & 0\\ 0 & 3.1875 & 0 & 0 & 0 & 0\\ 0 & 0 & 4.59 & 0 & 0 & 0\\ 0 & 0 & 0 & 2.04 & 0 & 0\\ 0 & 0 & 0 & 0 & 3.1875 & 0\\ 0 & 0 & 0 & 0 & 0 & 4.59\\ \end{bmatrix} \]

Simulando datos desde el modelo en R

\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]

# Matriz de covarianzas únicas poblacional
(Theta <- diag(c(2.04, 3.1875, 4.59, 2.04, 3.1875, 4.59), 6, 6))

     [,1]   [,2] [,3] [,4]   [,5] [,6]
[1,] 2.04 0.0000 0.00 0.00 0.0000 0.00
[2,] 0.00 3.1875 0.00 0.00 0.0000 0.00
[3,] 0.00 0.0000 4.59 0.00 0.0000 0.00
[4,] 0.00 0.0000 0.00 2.04 0.0000 0.00
[5,] 0.00 0.0000 0.00 0.00 3.1875 0.00
[6,] 0.00 0.0000 0.00 0.00 0.0000 4.59

Simulando datos desde el modelo en R

\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]

• ¡Ahora podemos calcular \(\color{#9a2515}{\boldsymbol\mu}\) y \(\color{#9a2515}{\mathbf{\Sigma}}\) tal y como indica el modelo!

# Parámetros del modelo multivariante
mu <- Nu    # Como alpha = 0, Lambda * alpha = 0.
Sigma = Lambda %*% Psi %*% t(Lambda) + Theta

• Finalmente, simulamos las respuesta de \(n=300\) participantes.

# Simulamos datos desde el modelo
library(MASS) # Tiene la distribución normal multivariante
set.seed(2025) # Simulación reproducible
Y <- mvrnorm(n = 300, mu = mu, Sigma = Sigma)
Y <- as.data.frame(Y) # Guardamos como data frame

Sintaxis del modelo

• La sintaxis de {blavaan} es idéntica a la de {lavaan}.
• Especificaremos dos modelos: uno con un único factor común y otro con dos factores relacionados (el verdadero).

# Sintaxis del modelo: un factor común
model.1f <- '
F1 =~ V1 + V2 + V3 + V4 + V5 + V6
'
# Sintaxis del modelo: dos factores relacionados
model.2f <- '
F1 =~ V1 + V2 + V3
F2 =~ V4 + V5 + V6
F1 ~~ F2
'

Prior predictive check

• ¿Qué priors está utilizando {blavaan} para los parámetros del modelo?

# By-default prior distributions
dpriors()

Prior predictive check

• {blavaan} tiene un atajo para realizar prior predictive checks.

# Ajusta un modelo "vacío" prisamp = TRUE!
prior_predictive <- bcfa(model.2f, data = Y, std.lv = TRUE, 
                         meanstructure = TRUE, test = "none", 
                         sample = 2e4, prisamp = TRUE, 
                         bcontrol = list(cores = 3)) 

• Acto seguido, podemos guardar la distribución a priori con la siguiente línea de código

# Guardamos la distribución a priori
prior.draws <- blavInspect(prior_predictive, "mcmc")

• Graficaremos las distribuciones con el paquete {bayesplot}

Prior predictive check: intersecciones

# Prior predictive distribution: intersecciones
mcmc_hist(x = prior.draws, regex_pars = "~1", alpha = .5)

Prior predictive check: pesos factoriales

# Prior predictive distribution: pesos factoriales
mcmc_hist(x = prior.draws, regex_pars = "=~", alpha = .5)

Prior predictive check: correlaciones latentes

# Prior predictive distribution: correlación latente
mcmc_hist(x = prior.draws, pars = "F1~~F2", alpha = .5)

Prior predictive check: varianzas únicas

# Prior predictive distribution: correlación latente
pars <- paste0("V", 1:6, "~~", "V", 1:6)
mcmc_hist(x = prior.draws, pars = pars,  alpha = .5)

Cambiando la distribución prior

• En {blavaan} no se puede cambiar la distribución a priori, sólo sus parámetros.

# Nueva distribución a priori
(new_priors <- dpriors(nu = "normal(0,15)", 
                       lambda = "normal(0,5)", 
                       rho = "beta(2,2)", 
                       psi = "gamma(1,.2)"))
# Muestreamos de la nueva distribución a priori
new_prior_pred <- bcfa(model.2f, data = Y, std.lv = TRUE, 
                       meanstructure = TRUE, test = "none", 
                       sample = 2e4, prisamp = TRUE, 
                       dp = new_priors, # Added new priors
                       bcontrol = list(cores = 3)) 

Cambiando la distribución prior

• Volvemos a extraer la distribución a priori

new.prior.draws <- blavInspect(new_prior_pred, "mcmc")

Prior predictive check: intersecciones

# Prior predictive distribution: intersecciones
mcmc_hist(x = new.prior.draws, regex_pars = "~1", alpha = .5)

Prior predictive check: pesos factoriales

# Prior predictive distribution: pesos factoriales
mcmc_hist(x = new.prior.draws, regex_pars = "=~", alpha = .5)

Prior predictive check: correlaciones latentes

# Prior predictive distribution: correlación latente
mcmc_hist(x = new.prior.draws, pars = "F1~~F2", alpha = .5)

Prior predictive check: varianzas únicas

# Prior predictive distribution: correlación latente
pars <- paste0("V", 1:6, "~~", "V", 1:6)
mcmc_hist(x = new.prior.draws, pars = pars,  alpha = .5)

Estimación del modelo

• Estimaremos los modelos de uno y dos factores con las nuevas priors.

• Utilizaremos 3 cadenas de Markov con 1500 iteraciones, descartando las primeras 500 iteraciones.

• bcontrol permite paralelizar las cadenas de Markov.

# Modelo de un factor
bcfa.1f.fit <- bcfa(model.1f, data = Y, burnin = 500, 
                    sample = 1000, meanstructure = TRUE, 
                    std.lv = TRUE, bcontrol = list(cores = 3))
# Modelo de dos factores
bcfa.2f.fit <- bcfa(model.2f, data = Y, burnin = 500, 
                    sample = 1000, meanstructure = TRUE, 
                    std.lv = TRUE, bcontrol = list(cores = 3))

Convergencia del modelo

• La validez de las inferencias depende de la convergencia del modelo.
• Potential Scale Reduction Factor (PSRF, \(\hat R\)) con valores inferiores a 1.05

# Potential Scale Reduction Factor (PSRF): 1 factor
which(blavInspect(bcfa.1f.fit, what = "rhat") > 1.05)

named integer(0)

# Potential Scale Reduction Factor (PSRF): 2 factores
which(blavInspect(bcfa.2f.fit, what = "rhat") > 1.05)

named integer(0)

• Esto sugiere que las cadenas de Markov están muestreando de la misma distribución posterior de manera estable.

Convergencia del modelo

• Podemos evaluar visualmente la convergencia de las cadenas de Markov a través de los gráficos de trayectoria (traceplots).
• Para ello, necesitamos extraer la distribución posterior de los parámetros.

# Distribución posterior de ambos modelos
posterior.1f <- blavInspect(bcfa.1f.fit, "mcmc")
posterior.2f <- blavInspect(bcfa.2f.fit, "mcmc")

• El siguiente código genera traceplots para los pesos factoriales

# Gráficos de trayectoria
mcmc_trace(x = posterior.1f, regex_pars = "F1=~")
mcmc_trace(x = posterior.2f, regex_pars = c("F1=~", "F2=~"))

Convergencia del modelo

Convergencia del modelo

Convergencia del modelo

• Podemos destacar una de las cadenas para evaluar su desempeño.

mcmc_trace_highlight(x = posterior.2f, regex_pars = "F1~~F2", 
                     highlight = 3)

Convergencia del modelo

• ¿Qué aspecto tiene este gráfico cuando la convergencia es mala?

Eficiencia

• Los valores de los parámetros de iteración a iteración suelen estar correlacionados.
• Cuanto menor sea su correlación, más informativo será el valor en cada muestra.
• Se evalúa a través del tamaño muestral efectivo.
  • Valor mínimo recomendable: \(100\times\) Nº de cadenas de Markov.

min(blavInspect(bcfa.1f.fit, "neff"))

[1] 947.9307

min(blavInspect(bcfa.2f.fit, "neff"))

[1] 2093.932

Ajuste del modelo

• Similares a los índices de ajuste frecuentistas, ¡pero con intervalos de credibilidad!
• Índices Bayesianos de ajuste global
  • Bayesian Root Mean Square Error of Approximation (BRMSEA)
  • Bayesian Unbiased Goodnes-of-fit Index (\(\Gamma\))
  • Bayesian Unbiased Goodnes-of-fit Index ajustado a muestras pequeñas (\(\Gamma_{adj}\))
  • Bayesian McDonald Index (BMc)
• Índices Bayesianos de ajuste incremental
  • Bayesian Comparative Fit Index (BCFI)
  • Bayesian Tucker-Lewis Index (BTLI)

Ajuste del modelo

# Sintaxis del modelo nulo
model.null <- paste0("V",1:6,"~~","V",1:6,collapse = " \n ")
# Ajustamos el modelo nulo
bcfa.null.fit <- bcfa(model.null, data = Y, burnin = 500, 
                      sample = 1000, meanstructure = TRUE, 
                      std.lv = TRUE, bcontrol = list(cores = 3))

Ajuste del modelo: un factor común

# Indices de ajuste: un factor
blav_fit_1f <- blavFitIndices(object = bcfa.1f.fit, 
                              baseline.model = bcfa.null.fit)
summary(blav_fit_1f, hpd = TRUE, prob = .95,
        central.tendency = c("mean","median","mode"))

               EAP Median   MAP    SD  lower upper
BRMSEA       0.313  0.313 0.312 0.007  0.301 0.327
BGammaHat    0.873  0.874 0.874 0.005  0.864 0.882
adjBGammaHat 0.228  0.230 0.235 0.030  0.169 0.282
BMc          0.804  0.805 0.806 0.008  0.789 0.818
BCFI         0.714  0.715 0.716 0.013  0.688 0.738
BTLI         0.052  0.056 0.059 0.044 -0.035 0.131
BNFI         0.714  0.715 0.716 0.013  0.688 0.737

Ajuste del modelo: un factor común

mcmc_hist(data.frame(blav_fit_1f@indices), alpha = 0.5,
          pars = c("BRMSEA", "BCFI", "BTLI", "BGammaHat"))

Ajuste del modelo: dos factores correlacionados

# Indices de ajuste: dos factores
blav_fit_2f <- blavFitIndices(object = bcfa.2f.fit, 
                              baseline.model = bcfa.null.fit)
summary(blav_fit_2f, hpd = TRUE, prob = .95,
        central.tendency = c("mean","median","mode"))

               EAP Median   MAP    SD lower upper
BRMSEA       0.046  0.050 0.057 0.026 0.000 0.085
BGammaHat    0.993  0.993 1.000 0.006 0.981 1.000
adjBGammaHat 0.975  0.978 0.998 0.021 0.936 1.000
BMc          0.989  0.990 0.999 0.009 0.972 1.000
BCFI         0.985  0.987 0.999 0.012 0.962 1.000
BTLI         0.974  0.976 0.983 0.024 0.924 1.017
BNFI         0.969  0.970 0.974 0.013 0.943 0.992

Ajuste del modelo: dos factores correlacionados

mcmc_hist(data.frame(blav_fit_2f@indices), alpha = 0.5,
          pars = c("BRMSEA", "BCFI", "BTLI", "BGammaHat"))

Comparación de modelos

• {blavaan} incluye dos métodos para comparar modelos:

1. Métodos basados en la capacidad predictiva del modelo:
   • Leave-One-Out (LOO) cross-validation.
   • Watanabee-Akaike Information Criteria (WAIC).
   • Ambos incluidos en el paquete {loo}.
   • Cuentan con pruebas de significacion.
2. Métodos basados en la capacidad explicativa del modelo:
   • Factor de Bayes (aproximación de Laplace-Metrópolis).
   • {blavaan} no recomienda su uso.

Comparación de modelos

# Comparación de modelos
blav_com_1f_2f <- blavCompare(object1 = bcfa.1f.fit, 
                              object2 = bcfa.2f.fit)
# Leave-One-Out: 
blav_com_1f_2f$loo[[1]] # modelo de un factor
blav_com_1f_2f$loo[[2]] # modelo de dos factores

Comparación de modelos

# Comparación de ELPDs
blav_com_1f_2f$diff_loo

       elpd_diff se_diff
model2   0.0       0.0  
model1 -63.9      12.2  

• Model2 es el modelo que hemos indicado previamente en object 2

# Comparación de modelos
blav_com_1f_2f <- blavCompare(object1 = bcfa.1f.fit, 
                              object2 = bcfa.2f.fit)

• Por tanto, el modelo de dos factores tiene mejor capacidad predictiva que el de un factor porque la diferencia \(\Delta\mbox{ELPD} = -63.9\) está fuera del rango \(\pm 2 \times 12.2\).

Estimaciones de los parámetros: medias posteriores

parameterEstimates(bcfa.2f.fit, standardized = TRUE, 
                   output = "text")

Latent Variables:
                   Estimate   Std.lv  Std.all
  F1 =~                                      
    V1                1.452    1.452    0.711
    V2                2.091    2.091    0.759
    V3                2.031    2.031    0.673
  F2 =~                                      
    V4                1.566    1.566    0.722
    V5                1.671    1.671    0.708
    V6                1.905    1.905    0.667

Covariances:
                   Estimate   Std.lv  Std.all
  F1 ~~                                      
    F2                0.469    0.469    0.469

Intercepts:
                   Estimate   Std.lv  Std.all
   .V1               -1.086   -1.086   -0.532
   .V2               -0.013   -0.013   -0.005
   .V3                1.029    1.029    0.341
   .V4               -0.918   -0.918   -0.423
   .V5                0.041    0.041    0.017
   .V6                0.986    0.986    0.345
    F1                0.000    0.000    0.000
    F2                0.000    0.000    0.000

Variances:
                   Estimate   Std.lv  Std.all
   .V1                2.060    2.060    0.494
   .V2                3.218    3.218    0.424
   .V3                4.974    4.974    0.547
   .V4                2.247    2.247    0.478
   .V5                2.779    2.779    0.499
   .V6                4.538    4.538    0.556
    F1                1.000    1.000    1.000
    F2                1.000    1.000    1.000

Distribución posterior: intersecciones

# Posterior distribution: histograms
mcmc_hist(x = posterior.2f, regex_pars = "~1", alpha = .5)

Distribución posterior: correlación latente

# Posterior distribution: density plot
mcmc_dens(x = posterior.2f, pars = "F1~~F2", alpha = .5)

Distribución posterior: pesos factoriales

# Posterior distribution: pairs plots
mcmc_pairs(posterior.2f, regex_pars = "F1=~", diag_fun = "hist", 
           off_diag_fun = "scatter", off_diag_args = list(alpha = 0.5))

Distribución posterior: varianzas únicas

# Posterior distribution: pairs plots
mcmc_intervals(posterior.2f, pars = paste0("V", 1:6, "~~", "V", 1:6))

Distribución posterior: pesos factoriales

# Posterior distribution: pairs plots
mcmc_areas(posterior.2f, regex_pars = c("F1=~", "F2=~"),
  prob = 0.8, prob_outer = 0.95, point_est = "median")

Posterior Predictive Model Checks (PPMC)

• Método para evaluar la capacidad predictiva de un modelo.
• La idea central es verificar si el model puede generar datos similares a los datos observados.
• ¿Podemos distinguir los datos observados de datos simulados desde el modelo?
• En {blavaan} podemos utilizar la función sampleData para generar datos desde el modelo.
• Tras simularlo, se comparan gráficamente los datos observados y simulados

Posterior Predictive Model Checks: \(y^{\text{rep}}\)

• Finalmente, podemos ver el parecido entre \(y\) con \(y^{\text{rep}}\)

# Simulamos 100 bases de datos
yrep <- sampleData(bcfa.2f.fit, nrep = 100, simplify = TRUE)
# Guardamos los resultados para cada item
posterior_plots_list <- vector(mode = "list", length = ncol(Y))
for(i in 1:ncol(Y)){
  posterior_plots_list[[i]] <- ppc_dens_overlay(
    y = Y[,i], yrep = t(sapply(yrep, function(x) x[,i]))
    ) + 
    labs(title = paste("Item", i))
}
# Gráfico comparando y e yrep
cowplot::plot_grid(plotlist = posterior_plots_list)

Posterior Predictive Model Checks: \(y^{\text{rep}}\)

Posterior Predictive Model Checks (PPMC)

• El PPMC también puede aplicarse a estadísticos específicos (SRMR, \(\chi^2\)).
• Durante la estimación, en cada iteración:
  1. Se extrae una muestra de los parámetros \(\theta^{(s)}\) de la distribución posterior (esto lo hacemos siempre).
  2. Se utilizan estos parámetros para generar datos simulados: \(y^{\text{rep}}\).
  3. Se calcula algún estadístico de interés en \(y^{\text{rep}}\) e \(y\).
  4. Se compara la distribución del estadístico entre \(y\) e \(y^{\text{rep}}\).
• En {blavaan} pueden realizarse con la función ppmc

Posterior Predictive Model Checks (PPMC)

• En el modelo bayesiano, en cada iteracion:
  • Estimamos \(\mathbf{\Lambda}^{(s)}\), \(\mathbf{\Theta}^{(s)}\) y \(\mathbf{\Psi}^{(s)}\).
  • Estimamos \(\mathbf{\Sigma}^{(s)} = \mathbf{\Lambda}^{(s)}\cdot\mathbf{\Psi}^{(s)}\cdot\mathbf{\Lambda}^{\prime(s)} + \mathbf{\Theta}^{(s)}\).
  • Siempre tenemos la misma matriz de covarianzas observada: \(\mbox{Cov}\left(Y\right)\).
  • Calculamos el estadístico \(\chi^{2(s)}\), que refleja la discrepancia entre \(\mbox{Cov}\left(Y\right)\) y \(\mathbf{\Sigma}^{(s)}\).
  • Simulamos \(y^{rep(s)} \sim \mathcal{MVN}\left(\boldsymbol{\mu}^{(s)}, \mathbf{\Sigma}^{(s)}\right)\)
  • Calculamos la media y covarianza “observada” de \(y^{rep(s)}\).
  • Calculamos el estadístico \(\chi^{2(s),rep}\), que refleja la discrepancia entre \(\mbox{Cov}\left(y^{rep(s)}\right)\) y \(\mathbf{\Sigma}^{(s)}\).

Posterior Predictive Model Checks (PPMC): \(\chi^2\) y SRMR

• La función ppmc puede utilizarse para calcular el estadístico \(\chi^2\) y SRMR en los datos observados y simulados en cada iteración.

# Posterior Predictive Model Checks
ppmc_fit_1f <- ppmc(object = bcfa.1f.fit, thin = 1,
                    fit.measures = c("srmr","chisq"))
ppmc_fit_2f <- ppmc(object = bcfa.2f.fit, thin = 1,
                    fit.measures = c("srmr","chisq"))

• Las funciones plot e hist comparan gráficamente los datos observados y simulados.

Posterior Predictive Model Checks (PPMC): \(\chi^2\)

plot(ppmc_fit_2f, element = "chisq") # Dos factores

Posterior Predictive Model Checks (PPMC): \(\chi^2\)

hist(ppmc_fit_2f, element = "chisq") # Dos factores

Posterior Predictive Model Checks (PPMC): \(\chi^2\)

plot(ppmc_fit_1f, element = "chisq") # Un factor

Posterior Predictive Model Checks (PPMC): \(\chi^2\)

hist(ppmc_fit_1f, element = "chisq")

Posterior Predictive Model Checks (PPMC): SRMR

plot(ppmc_fit_2f, element = "srmr") # Dos factores

Posterior Predictive Model Checks (PPMC): SRMR

hist(ppmc_fit_2f, element = "srmr") # Dos factores

Posterior Predictive Model Checks (PPMC): SRMR

plot(ppmc_fit_1f, element = "srmr") # Un factor

Posterior Predictive Model Checks (PPMC): SRMR

hist(ppmc_fit_1f, element = "srmr") # Un factor

PPMC y Fiabilidad: \(\omega\)

• Podemos utilizar ppmc con cualquier función que pueda utilizarse para un modelo de {lavaan}.
• Podemos utilizar compRelSEM de {SemTools} para obtener la distribución posterior de \(\omega\).

# Función para obtener la fiabilidad en cada iteración
bcfa.reliability <- function(fit){ semTools::compRelSEM(fit) }
# Fiabilidad en el modelo unidimensional
ppmc_omega_1f <- ppmc(object = blavaan.1f.fit, 
                      discFUN = bcfa.reliability)
# Fiabilidad en el modelo de dos factores
ppmc_omega_2f <- ppmc(object = blavaan.2f.fit, 
                      discFUN = bcfa.reliability)

PPMC y Fiabilidad: \(\omega\)

# Distribución posterior de la fiabilidad: un factor
hist(unlist(posterior_omega_1f@obsDist$discFUN1))

PPMC y Fiabilidad: \(\omega\)

Stan: state of the art

Stan: state of the art

• En R, podemos usar {Rstan} o {cmdstanr}.
• Es muy recomendable utilizar {cmdstanr}:
  1. La elección de los creadores de Stan.
  2. Sólo {cmdstanr} utiliza la última versión de Stan.
  3. Evita saturaciones en la memoria RAM.
  4. Más rápido y eficiente en la estimación.
  5. Conectado con otros paquetes: {posterior}, {loo}, {bayesplot}…
• De hecho, es más fácil instalar Stan usando {cmdstanr}.

Programación en Stan.

• Puede hacerse en Rstudio.
• Cuenta con siete posibles bloques:
  1. Functions: espacio para crear funciones específicas.
  2. Data: se introducen los datos desde R como una lista.
  3. Transformed data: transformaciones sobre los datos introducidos.
  4. Parameters: parámetros que muestrearemos de la posterior.
  5. Transformed parameters: transformaciones sobre los parámetros muestreados.
  6. Model: priors y log-verosimilitud.
  7. Generated quantities: parámetros transformados post-estimación.

Programación en Stan.

// Espacio para funciones
functions {
  
}
// Espacio para los datos
data {
  
}
// Espacio para transformar los datos
transformed data {

}
// Espacio para los parámetros
parameters {

}
// Espacio para los parámetros transformados
transformed parameters {

}
// Espacio para el modelo
model {

}
// Espacio para las cantidades generadas
generated quantities{

}

Programación en Stan: Data

• Programaremos el código para un modelo factorial confirmatorio.
• Los datos necesarios son:
  1. Número de participantes: \(N\)
  2. Número de ítems: \(J\)
  3. Número de latentes: \(M\)
  4. Matriz con los datos observados: \(Y\)
  5. Matriz que indique qué item satura en qué factor.
• Esta es mi forma de hacerlo, pero hay infinitas versiones alternativas!

Programación en Stan: Data

data {
  int<lower=0> N;             // Number of observations
  int<lower=0> J;             // Number of items
  int<lower=0> M;             // Number of latent factors
  matrix[N, J] Y;             // Observed data
  matrix[J, M] L_ind;         // Lambda index matrix
}

Programación en Stan: Parámetros

• Vamos a especificar que el modelo muestree
  1. El vector de medias del modelo: \(\boldsymbol\mu\).
  2. La matriz de pesos factoriales: \(\mathbf{\Lambda}\).
  3. El vector de varianzas únicas: \(\boldsymbol\theta\).
  4. La matriz triangular inferior de la descomposición de Cholesky para las correlaciones latentes: \(L_\boldsymbol\Psi\)
     • Garantiza que la matriz de correlaciones es válida
     • Recuperamos \(\boldsymbol\Psi\) como \(\boldsymbol\Psi=L_\boldsymbol\Psi\cdot L_\boldsymbol\Psi^\prime\)

Programación en Stan: Parámetros

parameters {
  vector[J] mu;
  matrix[J, M] L_unc;
  vector<lower=0>[J] Theta;
  cholesky_factor_corr[M] Psi_chol;
}

Programación en Stan: Parámetros transformados

1. Indicar qué peso va con cada factor latente.
2. Construir \(\mathbf{\Psi}\) como \(L_\mathbf{\Psi}\cdot L_\mathbf{\Psi}^\prime\).
3. Construir \(\mathbf{\Sigma}\) como \(\mathbf{\Sigma} = \mathbf{\Lambda}\cdot\mathbf{\Psi}\cdot\mathbf{\Lambda}^\prime+\mathbf{\Theta}\) .

transformed parameters{
  // Constrain Lambda matrix
  matrix[J, M] L_mat = L_unc .* L_ind;
  // Model-implied latent-variable correlation matrix
  corr_matrix[M] Psi_b = Psi_chol + Psi_chol';
  // Model-implied covariance matrix
  cov_matrix[J] Sigma;
  Sigma = L_mat * Psi_b * L_mat' + diag_matrix(Theta);
}

Programación en Stan: Modelo

• Sólo debemos poner priors en los parámetros muestreados.
• En este ejemplo usaremos priors poco informativas.

• Finalmente, especificaremos la función de verosimilitud del modelo.

Programación en Stan: Modelo

model {
  // Model priors
  mu ~ normal(0, 10);              // Mu prior
  to_vector(L_unc) ~ normal(0, 5); // Lambda prior
  Theta ~ gamma(1, 0.5);           // Theta prior
  Psi_chol ~ lkj_corr_cholesky(1); // Psi prior
  
  // Model log-likelihood
  for(i in 1:N){
    Y[i,] ~ multi_normal(mu, Sigma);
  }
}

Programación en Stan: Cantidades generadas

• Este apartado permite realizar cálculos tras ajustar el modelo.
  • Otras transformaciones de parámetros.
  • Simulación de datos (\(y^{\text{rep}}\))
  • Posterior Predictive Model Checks
  • Cálculo de log-verosimilitudes
• Nosotros guardaremos la log-verosimilitud para comparar modelos.

generated quantities {
  vector[N] log_lik;          // Model log-likelihood
  for(i in 1:N){
    log_lik[i] = multi_normal_lpdf(Y[i,] | mu, Sigma);
  }
}

Programación en Stan: Cantidades generadas

• El modelo estima los parámetros \(\boldsymbol\mu\) y \(\mathbf\Sigma\) en función de los datos, \(y\).
• Sin embargo, podemos obtener la misma \(\mathbf\Sigma\) con distintas \(\mathbf\Lambda\):
  • Obtendríamos la misma \(\mathbf\Sigma\) usando \(\mathbf\Lambda\) o \(-\mathbf\Lambda\).
• {blavaan} soluciona este problema de signos en las cantidades generadas.
• El código de este modelo incluye esta corrección.
  • No lo mostraremos en detalle, basta con saber que se aplica para que el modelo esté identificado.

Usando {cmdstanr}

• El primer paso es compilar nuestro modelo programado en Stan

# Compilamos el modelo
BCFA <- cmdstan_model(stan_file = "CFA_marginal.stan")

• Una vez compilado, {cmdstanr} creará un programa ejecutable en vuestro ordenador.
• Al mismo tiempo, tendremos acceso a este programa con el objeto donde hayamos guardado el modelo en R.
• El objeto BCFA es una especie de lista con todas las funciones que pueden utilizarse dentro.

Usando {cmdstanr}

• Por ejemplo, BCFA$sample() es la función para estimar el modelo bayesiano.
• Para introducir los datos en el modelo, utilizaremos una lista.
• Los objetos de la lista deben coincidir en nombre con los del apartado de Datos en Stan

sdata.1f <- list(
  N = 300,
  J = 6,
  M = 1,
  Y = Y,
  L_ind = matrix(rep(1, 6),ncol=1)
)

Usando {cmdstanr}

• BCFA$sample() cuenta con muchos argumentos, pero no son similares a los de {blavaan}

BCFA_1f <- BCFA$sample(
  data = sdata.1f,       # Stan data
  chains = 4,            # Number of chains
  parallel_chains = 4,   # Number of parallel chains
  iter_warmup = 500,     # Adaptation iterations
  iter_sampling = 1500,  # Sampled iterations
  refresh = 500,         # Progress bar at 500 iterations
  init = 0)              # All starting values = 0

• BCFA-1f será una lista que contenga los resultados y funciones para almacenarlos y resumirlos.

Usando {cmdstanr}

• La función BCFA_1f$summary() genera tablas resumen de los parámetros que indiquemos.

BCFA_1f$summary("Lambda") # Estimated factor loadings

Usando {cmdstanr}

• Estimaremos también el modelo de dos factores.

sdata.2f <- list(
  N = 300,
  J = 6,
  M = 2,
  Y = Y,
  L_ind = ifelse(Lambda == 0, 0, 1)
)

Usando {cmdstanr}

BCFA_2f <- BCFA$sample(
  data = sdata.2f,       # Stan data
  chains = 4,            # Number of chains
  parallel_chains = 4,   # Number of parallel chains
  iter_warmup = 500,     # Adaptation iterations
  iter_sampling = 1500,  # Sampled iterations
  refresh = 500,         # Progress bar
  seed = 2025,           # Reproducible results
  init = 0)              # All starting values = 0

Usando {cmdstanr}

BCFA_2f$summary("Lambda") # Estimated factor loadings

Comparando modelos

• Hemos guardado la log-verosimilitud de ambos modelos.
• Podemos acceder a todos sus valores usando $draws()
• Los utilizaremos para comparar los modelos con {loo}

loo_1f  <- loo(x = BCFA_1f$draws("log_lik"))
loo_2f  <- loo(x = BCFA_2f$draws("log_lik"))
loo_compare(loo_1f, loo_2f)

       elpd_diff se_diff
model2   0.0       0.0  
model1 -55.2      12.0  

Comparación de modelos

Comparación de modelos: el Factor de Bayes

• La verosimilitud marginal es la probabilidad de los datos observados con el modelo que hemos ajustado.

• El Factor de Bayes es la ratio entre la verosimilitud marginal de dos modelos \(\mathcal{M}\).

Comparación de modelos: el Factor de Bayes

\[ \mbox{BF}_{12}=\frac{\Pr(D\mid\mathcal{M}_1)}{\Pr(D\mid\mathcal{M}_2)} \]

• El Factor de Bayes compara la probabilidad de los datos observados en cada modelo.
• Permite responder a la pregunta: ¿Con qué modelo es más probable que se hayan generado los datos observados?.
• \(\mbox{BF}_{12}=3\) significa que los datos observados son tres veces más probables con el modelo \(\mathcal{M}_1\) que con \(\mathcal{M}_2\).
• Su cálculo exacto no está disponible en modelos complejos, pero puede aproximarse con métodos numéricos.

Limitaciones del factor de Bayes.

• Las críticas más relevantes contra el Factor de Bayes son:
  1. Las aproximaciones pueden ser inestables y sesgadas en modelos complejos.
  2. Es muy sensible a las distribuciones a priori.
  3. Puede requerir muchas extracciones de la posterior.
• Existen recomendaciones para:
  1. Evaluar la inestabilidad y posible sesgo del Factor de Bayes.
  2. Evaluar el impacto de la distribución a priori.

Comparación de modelos: Leave-One-Out

• El método Leave-One-Out (LOO) evalúa la capacidad de un modelo para predecir datos que no fueron utilizados en su ajuste.
• Permite responder a la pregunta: ¿Cómo de buenas son las predicciones del modelo fuera de la muestra original?
• El procedimiento es sencillo:
  1. Quitamos una observación de la base de datos.
  2. Ajustamos el modelo con los datos restantes.
  3. Predecimos la observación que hemos quitado al principio.
  4. Calculamos la log-densidad predictiva
  5. Repetimos esto para todas las observaciones.
  6. Sumamos las log-densidades predictivas.

Comparación de modelos: Leave-One-Out

• Ejemplo: regresión lineal simple (\(Y_i = \beta_0+\beta_1\cdot X_i + \varepsilon_i\))

for (i in 1:n) {
  # Ajuste sin la i-ésima observación
  fit_loo <- lm(y[-i] ~ x[-i])
  # Predicción para la observación i
  y_pred_i <- predict(fit_loo, newdata = data.frame(x = x[i]))
  # Error de estimación sin la i-ésima obs
  sigma_i <- summary(fit_loo)$sigma
  # Cálculo del log-densidad predictiva de y[i]
  log_pdens[i] <- dnorm(y[i], mean = y_pred_i, sd = sigma_i, 
                        log = TRUE)
}
# ELPD: suma de las log-densidades
ELPD_loo <- sum(log_pdens)

Comparación de modelos: Leave-One-Out

• ELPD son las siglas de Expected log pointwise predictive density, y su valor refleja la calidad predictiva del modelo.

• La comparación entre ELPDs permite saber qué modelo predice mejor los datos fuera de la muestra.

• El paquete {loo} estima \(\Delta\mbox{ELPD}_{\mathcal{M}_1,\,\mathcal{M}_2}\) y su error típico.
• Si \(\Delta\mbox{ELPD}_{\mathcal{M}_1,\,\mathcal{M}_2}\) es mayor que \(\pm 2\, \times\) el error típico, la diferencia predictiva entre ambos modelos es significativa.

Limitaciones de Leave-One-Out

• Las críticas más relevantes son:
  • Computacionalmente costoso, intratable en bases de datos grandes.
  • Sigue siendo sensible a las priors, aunque menos que el \(BF\).
  • Asumen que los datos son intercambiables (i.e., no hay un efecto del tiempo).
• Hay alternativas para algunos de los problenas:
  • PSIS-LOO: versión adaptada a bases de datos grandes.
  • Leave-future-one-out para tener en cuenta la secuencia temporal.

Conectando ambos mundos: Pseudo-factor de Bayes

• El Pseudo-factor de Bayes (\(\mbox{PsBF}\)) compara la capacidad predictiva de dos modelos de forma similar al Factor de Bayes.
• Permite responder a la pregunta: ¿Qué modelo predice mejor los datos observados?

• \(\mbox{PsBF}_{12}=3\) significa que las predicciones con el modelo \(\mathcal{M}_1\) son tres veces mejores que las predicciones con el modelo \(\mathcal{M}_2\).

Comparación de modelos

• {blavaan} incluye una estimación de la verosimilitud marginal y del factor de bayes basada en la aproximación de Laplace-Metrópolis.

# Comparación de ELPDs
blav_com_1f_2f$bf

         bf        mll1        mll2 
  -62.32224 -4102.84295 -4040.52072 

• Como está en escala logarítmica, debemos exponenciarlo.

c(BF12 = exp(blav_com_1f_2f$bf[1]), # BF12
  BF21 = 1/exp(blav_com_1f_2f$bf[1])) # BF21

     BF12.bf      BF21.bf 
8.586119e-28 1.164671e+27 

Comparación de modelos

• Los datos observados son \(1.16\times 10^{27}\) veces más probables con el modelo de dos factores que con el modelo de un factor.

• Finalmente, el Pseudo-Factor de Bayes se obtiene exponenciando \(\Delta\mbox{ELPD}\).

# Incremento el ELPD
ELPD_diff <- blav_com_1f_2f$diff_loo[2,"elpd_diff"]
c(PsBF12 = exp(ELPD_diff), PsBF21 = 1/exp(ELPD_diff))

      PsBF12       PsBF21 
1.761748e-28 5.676181e+27 

• La capacidad predictive del modelo de dos factores es \(5.67\times 10^{27}\) veces mejor que la del modelo de un factor.